Magnetohydrodynamics in full general relativity: Formulation and tests 
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A new implementation for magnetohydrodynamics (MHD) simulations in full general relativity 
(involving dynamical spacetimes) is presented. In our implementation, Einstein's evolution equations 
are evolved by a BSSN formalism, MHD equations by a high-resolution central scheme, and induction 
equation by a constraint transport method. We perform numerical simulations for standard test 
problems in relativistic MHD, including special relativistic magnetized shocks, general relativistic 
magnetized Bondi flow in stationary spacetime, and a longterm evolution for self-gravitating system 
composed of a neutron star and a magnetized disk in full general relativity. In the final test, we 
illustrate that our implementation can follow winding-up of the magnetic field lines of magnetized 
and differentially rotating accretion disks around a compact object until saturation, after which 
magnetically driven wind and angular momentum transport inside the disk turn on. 
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I. INTRODUCTION 



Hydrodynamics simulation in general relativity is probably the best theoretical approach for investigating dynamical 
phenomena in relativistic astrophysics such as stellar core collapse to a neutron star and a black hole, and the 
merger of binary neutron stars. In the past several years, this field has been extensively developed (e.g., [1-7]) and, 
as a result, now it is feasible to perform accurate simulations of such general relativistic phenomena for yielding 
scientific results (e.g., [6-9] for our latest results). For example, with the current implementation, radiation reaction 
of gravitational waves in the merger of binary neutron stars can be taken into account within ~ 1% error in an 
, appropriate computational setting [6,7]. This fact illustrates that the numerical relativity is a robust approach for 
■ detailed theoretical study of astrophysical phenomena and gravitational waves emitted. 

However, so far, most of the scientific simulations in full general relativity have been performed without taking into 
account detailed effects except for general relativistic gravity and pure hydrodynamics. For example, simplified ideal 
equations of state have been adopted instead of realistic ones (but see [7]). Also, the effect of magnetic fields has 
been neglected although it could often play an important role in the astrophysical phenomena (but see [10]). In the 
next stage of numerical relativity, it is necessary to incorporate these effects for more realistic simulations. As a step 
toward a more realistic simulation, we have incorporated an implementation for ideal magnetohydrodynamics (MHD) 
r* ' equations in fully general relativistic manner. In this paper, we describe our approach for these equations and then 
present numerical results for test problems computed by the new implementation. 

Magnetic fields indeed play an important role in determining the evolution of a number of relativistic objects. In 
the astrophysical context, the plasma is usually highly conducting, and hence, the magnetic fields are frozen in the 
matter. This implies that a small seed field can wind up and grow in the complex motion of the matter, resulting 
in a significant effect in the dynamics of the matter such as magnetically driven wind or jet and angular momentum 
redistribution. Specifically, in the context of the general relativistic astrophysics, the magnetic fields will play a role in 
the following phenomena and objects: Stellar core collapse of magnetized massive stars to a protoneutron star [11] or a 
black hole, stability of accretion disks (which are either non-self-gravitating or self-gravitating) around black holes and 
neutron stars, magnetic braking of differentially rotating neutron stars [10] which are formed after merger of binary 
neutron stars [6,7] and stellar core collapse [14-16,8,9], and magnetically induced jet around the compact objects 
(e.g., [17]). To clarify these phenomena, fully general relativistic MHD (GRMHD) simulation (involving dynamical 
spacetimes) is probably the best theoretical approach. 

In the past decade, numerical implementations for GRMHD simulation in the fixed gravitational field have been 
extensively developed (e.g., [18,17,19-24]). In particular, it is worth to mention that Refs. [19-23] have recently 
presented implementations for which detailed tests have been carried out for confirmation of the reliability of their 
computation, in contrast with the attitude in an early work [17]. They are applied for simulating magnetorotational 
instability (MRI) of accretion disks and subsequently induced winds and jets around black holes and neutron stars. 
On the other hand, little effort has been paid to numerical implementations of fully GRMHD (in the dynamical 
gravitational field). About 30 years ago, Wilson performed a simulation for collapse of a magnetized star in the 
presence of poloidal magnetic fields in general relativity. However, he assumes that the three-metric is conformally 
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flat [25], and hence, the simulation is not fully general rclativistic, although recent works have indicated that the 
conformally flat approximation works well in the axisymmctric collapse (e.g., compare results among [15], [26], and 
[27]). The first fully GRMHD simulation for stellar collapse was performed by Nakamura about 20 years ago [28]. He 
simulated collapse of nonrotating stars with poloidal magnetic fields to investigate the criteria for formation of black 
holes and naked singularities. Very recently, Duez et al. have presented a new implementation capable of evolution 
for the Einstein-Maxwell-MHD equations for general cases [10]. They report successful results for test simulations. 
Valencia group has also developed a GRMHD implementation very recently [29] . 

In this paper, we present our new implementation for fully GRMHD which is similar to but in part different from 
that in [10] *. As a first step toward scientific simulations, we have performed simulations in standard test problems 
including special relativistic magnetized shocks, general relativistic Bondi flow in stationary spacetime, and long term 
evolution of fully general relativistic stars with magnetic fields. We here report the successful results for these test 
problems. 

Before proceeding, we emphasize that it is important to develop new GRMHD implementations. In the presence of 
magnetic fields, matter motion often becomes turbulence-like due to MRIs in which a small scale structure often grows 
most effectively [30]. Furthermore in the presence of general relativistic self-gravity which has a nonlinear nature, 
the matter motion may be even complicated. Perhaps, the outputs from the simulations will contain results which 
have not been well understood yet, and thus, are rich in new physics. Obviously high accuracy is required for such 
frontier simulation to confirm novel numerical results. However, because of the restriction of computational resources, 
it is often very difficult to get a well-resolved and completely convergent numerical result in fully general relativistic 
simulation. In such case, comparison among various results obtained by different numerical implementations is crucial 
for checking the reliability of the numerical results. From this point of view, it is important to develop several numerical 
implementations in the community of numerical relativity. By comparing several results computed by different 
implementations, reliability of the numerical results will be improved each other. Our implementation presented 
here will be useful not only for finding new physics but also for checking numerical results by other implementations 
such as that very recently presented in [10,29]. 

In Sec. II, we present formulations for Einstein, Maxwell, and GRMHD equations. In Sec. Ill, numerical methods 
for solving GRMHD equations are described. In Sec. IV, methods for a solution of initial value problem in general 
relativity is presented. In Sees. V and VI, numerical results for special and general relativistic test simulations are 
shown. In the final subsection of Sec. VI, we illustrate that our implementation can follow growth of magnetic fields of 
accretion disks in fully general relativistic simulation. Sec. VII is devoted to a summary and a discussion. Throughout 
this paper, we adopt the geometrical units in which G = c = 1 where G and c are the gravitational constant and the 
speed of light. Latin and Greek indices denote spatial components and spacetime components, respectively, r^„ and 
Sij(= denote the flat spacetime metric and the Kronecker delta, respectively. 



II. BASIC EQUATIONS 
A. Definition of variables 

Basic equations consist of the Einstein equations, general relativistic hydrodynamic equations, and Maxwell equa- 
tions. In this subsection, we define the variables used in these equations. The fundamental variables for geometry are 
a: lapse function, [3 k : shift vector, 7^ : metric in three-dimensional spatial hypersurface, and Kif extrinsic curvature. 
The spacetime metric g M „ is written as 

lnv n^^ui (1) 
where is a unit normal to a spacelike spatial hypersurface S and is written as 

«" = (-,--), or n„ = (-a,0). (2) 

In the BSSN formalism [35], one defines 7 = rje 12 ^ — det(7y): determinant of 7^ , 7^ = e^^-jij: conformal three- 
metric, K — K k k : trace of the extrinsic curvature, and = e^ 4 ^(Kij — Kjij/3): a tracefree part of the extrinsic 



*For instance, our formulation for Einstein's evolution equations, gauge conditions, and our numerical scheme for GRMHD 
equations are different from those in [10] as mentioned in Sees. II and III. 



2 



curvature. Here, rj denotes the determinant of flat metric; in the Cartesian coordinates, r\ = 1, and in the cylindrical 
coordinates (zu, tp,z), r\ — w 2 . In the following, V M , Di, and Di denote the covariant derivatives with respect to <? M „, 
-fij, and jij, respectively. A and A denote the Laplacians with respect to 7^ and 7^ . Rij and Rij denote the Ricci 
tensors with respect to 7^ and jij , respectively. 

The fundamental variables in hydrodynamics are p: rest-mass density, e : specific internal energy, P : pressure, and 
: four velocity. From these variables, we define the following variables which often appear in the basic equations: 



p* 


= pwe 6 ^, 






_ dx l 


u l 




~~dt ~ 


u* 


h 




P 


= 1 + £ + 








P 


w 


= auK 





(3) 

Uj 

u' 



-F+J^, (4) 



(5) 
(6) 

Here, p» is a weighted baryon rest mass density from which the conserved baryon rest mass can be computed as 

M* = / prf^tPx. (7) 



The fundamental variable in the ideal MHD is only b^: magnetic field. The electric field in the comoving frame 
F^ v u v is assumed to be zero, and electric current j M is not explicitly necessary for evolving the field variables. Using 
the electromagnetic tensor F^ v , is defined by [31] 

b„ = -\e^ va pu v F a P, (8) 

where e^ap is the Levi-Civita tensor with etm = \/—g and e* 123 = —1/y/^g. Equation (8) implies 

= 0. (9) 

Using Eq. (8), F^ v in the ideal MHD is written as 
and thus, it satisfies the ideal MHD condition 

F^u" = 0. (11) 

The dual tensor of F^ is defined by 

F% = \^™pF at) = b^u v - b uUlx . (12) 

For rewriting the induction equation for the magnetic fields into a simple form (see Sec. II D), we define the three- 
magnetic field as 

B l = -e^F*^ = ^{wb 1 - atfu 1 ). (13) 
Here, we note that B* = (i.e., B^n^ = 0), and thus, Bi — JijB^. Equations (13) and (9) lead to 

Using the hydrodynamic and electromagnetic variables, energy-momentum tensor is written as 

T lu ,=T%T l + T™ (15) 
rJJ uld and T^ 1 denote the fluid and electromagnetic parts defined by 

^ Uid = (P + P£ + + P 9llv = phu^u v + P 9llv , (16) 



rpEM _ rp j? a _ 
± fj,i> — r liar v 



\g^F a0 F a ^ = Q 5aw + u^b 2 - b^K, (17) 
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where 

Thus, is written as 

T M „ = (ph + b 2 )u^u v + (p + \b^)g^ - b^K. (19) 

For the following, we define magnetic pressure and total pressure as P ma g = b 2 /2 and Ptot = P + b 2 /2, respectively. 
The (3+1) decomposition of T Mi , is 

p H = T^n" = (ph + b 2 )w 2 - P tot - (ab*) 2 , (20) 

Ji = -T^ 1 v l = (ph + b 2 )w Ul - ab%, (21) 

Sij = T^cfj = (P h + b2 ) u i u J + p totlij - hbj. (22) 
Using these, the energy-momentum tensor is rewritten in the form 

T^v = Pun^ + Ji-f\n v + Ja^n^ + Sirf^^. (23) 

This form of the energy-momentum tensor is useful for deriving the basic equations for GRMHD presented in Sec. II 
C. For the following, we define 

So = e 6< W (24) 
S t = e^Ji. (25) 

These variables together with p* and B l are evolved explicitly in the numerical simulation of the ideal MHD (see Sec. 
II C). 



B. Einstein's equation 

Our formulation for Einstein's equations is the same as in [6] in three spatial dimensions and in [34] in axial 
symmetry. Here, we briefly review the basic equations in our formulation. Einstein's equations are split into constraint 
and evolution equations. The Hamiltonian and momentum constraint equations are written as 

R k k - AijAv + ^K 2 = 167TPH, (26) 

DiA^-^DjK = 87vJ j: (27) 

or, equivalently 

= |p fc fe _ 2 W 5 - ^ (A^ 2 -K 2 ) , (28) 

AOA 6 ^,) - \^DjK = SttJ^ 6 , (29) 

where tp = e^. These constraint equations are solved to set initial conditions. A method in the case of GRMHD is 
presented in Sec. IV. 

In the following of this subsection, we assume that Einstein's equations are solved in the Cartesian coordinates 
(x,y,z) for simplicity. Although we apply the implementation described here to axisymmetric issues as well as 
nonaxisymmetric ones, this causes no problem since Einstein's equations in axial symmetry can be solved using the 
so-called Cartoon method in which an axisymmetric boundary condition is appropriately imposed in the Cartesian 
coordinates [32-34]: In the Cartoon method, the field equations are solved only in the y = plane, and grid points 
of y = ±Ax (Ax denotes the grid spacing in the uniform grid) are used for imposing the axisymmetric boundary 
conditions. 
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We solve Einstein's evolution equations in our latest BSSN formalism [35,6]. In this formalism, a set of variables 
(■jij, <j>, Aij, K, Fi) are evolved. Here, we adopt an auxiliary variable Fi = S^dijij that is the one originally proposed 
and different from the variable adopted in [10] in which is used. Evolution equations for 7^ , <f), , and K are 



(dt-/3 l di)A tJ =e 



■( 



R 



k 

„ -e^R k K ) - ( DiD ja - -e^%Aa 



(30) 



+a(KA ij - 2A ik A k ) + p k ^A k] + f3 k ^A ki - -p\ k A t3 
{d t ~p l d l )$= l -(-aK + p k M ) 1 



AijAv + -K 2 



Aa + Aira( PB , + S k k ). 



(d t -l3 l di)K = u 

For a solution of <j>, the following conservative form may be adopted [6] : 

d t e H -d^e^) = -aKe 6<t> . 
For computation of Ri 3 in the evolution equation of A i3 , we decompose 



Ri 



where 



R} 



-2DiDj<f> - 2%A(j) + ADi4>D j( j) - A%D k <t>D k 



Rij - 2 



S kl (-hijM + hik,ij + hjk,u) + 2d k (f kl Ti ti3 ) - 2f^f^ 



(31) 
(32) 
(33) 

(34) 

(35) 

(36) 
(37) 



In Eq. (37), we split 7^ and 7^ as 8i 3 + hi 3 and (5 y + / y , respectively, f k 3 is the Christoffel symbol with respect to 
jij, and Tk,ij — lki^\ 3 - Because of the definition det(7y) = 1 (in the Cartesian coordinates), we use T ki — 0. 

In addition to a flat Laplacian of h i3 , R i3 involves terms linear in hi 3 as 8 kl h ik ^ 3 + 8 kl h 3k _u- To perform numerical 
simulation stably, we replace these terms by Fij + F 3 ^. This is the most important part in the BSSN formalism, 
pointed out originally by Nakamura [28] . The evolution equation of Fi is derived by substituting Eq. (30) into the 
momentum constraint as 



{dt - P l di)Fi = -WnaJi + 2a{f kj A ikd + f k \ 3 A lk - -A?% u + §<P, k A\ - -K^ 
+ V k \-2 aM A l3 +{3\ k h l3 , + (7^. + 7jl f3 l ti - paPj) J. 



(38) 



We also have two additional notes for handling the evolution equation of A i3 . One is on the method for evaluation 
of R k for which there are two options, use of the Hamiltonian constraint and direct calculation by 



Rij7 ij =e- 4 *(R k + Rl j fi). 



(39) 



We always adopt the latter one since with this, the conservation of the relation A^j^ = is much better preserved. 
The other is on the handling of a term of -y 4 - 7 5 kl h i3 ^ k i which appears in R k . This term is written by 



f j 6 kl h ijM = -5 kl h lhk f\ v 



(40) 



where we use det(7,j) = 1 (in the Cartesian coordinates). 

As the time slicing condition, an approximate maximal slice condition K m is adopted following previous papers 
(e.g., [36]). As the spatial gauge condition, we adopt a hyperbolic gauge condition as in [37,6]. Successful numerical 
results for merger of binary neutron stars and stellar core collapse in these gauge conditions are presented in [6,7,26,8]. 
We note that these are also different from those in [10]. 



5 



C. GRMHD equations 



Hydrodynamic equations are composed of 



V^pu") = 0, 
7 /V M Tt = 0, 
n"V M Tt = 0. 



(41) 
(42) 
(43) 



The first, second, and third equations are the continuity, Euler, and energy equations, respectively. In the following, the 
equations are described for general coordinate systems since the hydrodynamic equations are solved in the cylindrical 
coordinates as well as in the Cartesian coordinates. 
The continuity equation (41) is immediately written to 



d t p* + -^pdi(p*y/rjv l ) = 0. 



Equations (42) and (43) are rewritten as 
Then, using Eq. (23), they are written to 



-T^dig^ = 0, 



g(pnn» + 7 " i J i )] = V7( a^Sy - J<D'a 



V7( -Padja + JidjF - ^S lk 8 ol lk 



where we use 



l n v djg ilv 



-2dj In a, 



V M n„ = -K^ v - n^Djj In a. 
The explicit forms of Eqs. (47) and (48) are 

d t S 3 + -U [Vv{sy + ae^P tot 6/ - ^B\B, + u 3 B k u k )}] 



= -Sodja + Skdjp 



k ' ae 6 * 



2S k d 3 



P tot d j In ^rj 



ae 2 *S ik drf k , 



Vv 

= lae°+KS k k + ae 2 *S t1 A 13 - S k D h 



a, 



where 



Si j Sij Ptot^fij- 



(44) 



(45) 
(46) 



(47) 
(48) 



(49) 
(50) 
(51) 
(52) 



(53) 



(54) 



(55) 



In the axisymmetric case, the equations for (p», Si, Sq) should be written in the cylindrical coordinates {vj,ip,z) 
when we adopt the Cartoon method for solving Einstein's evolution equations [32-34]. On the other hand, in the 
standard Cartoon method, Einstein's equations are solved in the y = plane for which x = zu, S w = S x , S v = xS y , 
and other similar relations hold for vector and tensor quantities. Taking into this fact, the hydrodynamic equations 
in axisymmetric spacetimes may be written using the Cartesian coordinates replacing (w,<p) by (x,y). Then, 
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dtp* + -d x {p*v x x) + d z (p*v z ) = 0, 



(56) 



d t S A + -d x 
x 



{s A v x + ae 6 *P tot 5 A x - 
S A v z + ae^Ptotd/ - 



a 



B x (b a + u a B 1 u^ 



+ 



S A a + S k d A f3 k + ae 6 * 
S„v y a 



L a 



B z (B A + u A B l u t 
2S k k d A 4> + ^S A X 



°A > 



-ae 2 *S lk d A f k 



-B v {B y + B l Ui u y ) 
^ + ±d x [x 2 { Sy v x - -^B x (B y + u y Wu) } 



+ d z [s y v z ^-^B z (B y + u y B^) 



= 0, 



dtS + \d x [x{s Q v x + e^P tot (v x + p x ) - -^B l u % B x } 
= \ae 6rt, KS k k + ae 2 ^S l3 A 11 - S k D k a, 



S v z + e^P tot (v z + (3 Z ) - -^B l Ul B> 



(57) 



(58) 



(59) 



where A denotes x or z, while i, j,k, ■ ■ ■ are x or y or z. 

After evolving p*, Si, and So together with B l (see next subsection for the equations), we have to determine the 
primitive variables such as p, e, Ui, and u* (or w = au*). For this procedure, we make an equation from the definition 

of Si as 



s 2 = P': 2 1 13 S 1 S 3 = (h+ B 2 w- 1 ^j 2 {w 2 - 1) - D 2 (hw)- 2 (B 2 + 2hw), 



where B 2 and D 2 are determined from the evolved variables (p*, Si,B l , <f>) as 



B 



B 2 



p*e 



and D 2 



{B l S t f 
pie** 



(60) 



(61) 



and for getting Eq. (60), we use the relation SiB 1 = p Sf hB % u i . Equation (60) is regarded as a function of h and w for 
given data sets of s 2 , B 2 , and D 2 . From the definition of So, we also make a function of h and w as 



^ = hw -JL + B 2 - -L(B 2 + D 2 h- 2 ). 
p* pw 2w z 



(62) 



Here, Pj p may be regarded as a function of h and w for a given data sets of p* and So. This is indeed the case for 
frequently used equations of state such as T-law equations of state P = (T — l)pe where T is an adiabatic constant and 
hybrid equations of state for which P is written in the form P(p, h) (e.g., [15,7]). Thus, Eqs. (60) and (62) constitute 
simultaneous equations for h and w for given values of p*, Si, So, B l , and geometric variables. The solutions for h 
and w are numerically computed by the Newton-Raphson method very simply. Typically, a convergent solution is 
obtained with four iterations according to our numerical experiments. 



D. Maxwell equations 

The Maxwell equations are 

V M F^ = -Anf, (63) 
V li F aP + V a F Pli + VpF lta =0. (64) 

In the ideal MHD, Eq. (63) is not necessary, and only Eq. (64) has to be solved. Using the dual tensor, Eq. (64) is 
rewritten to 

v m f; m = 0- (65) 

This immediately leads to 
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d k (^ 2 B k )=0, 



d t B k = 



1/2 



n 



V 1/2 {B i v k -B k v l ) 



(66) 
(67) 



Equation (66) is the no-monopoles constraint, and Eq. (67) is the induction equation. The constraint equation (66) 
is solved in giving initial conditions, and the induction equation is solved for the evolution. 
In the axisymmetric case, these equations in the y = plane are written as 



-d x (xB x ) + d z B z = 0, 
x 



d t B x = -d z (B x v z - B z v x ), 



(68) 
(69) 
(70) 

d t B y = d x (B x v y - B y v x ) + d z (B z v y - B y v z ). (71) 
Equations (69)-(71) together with Eqs. (56)-(59) constitute basic equations for ideal MHD in the axisymmetric case. 



d t B z = -a 



x(B x v z - B z v x ) 



E. Definition of global quantities 

In numerical simulations for self-gravitating system, in addition to the total baryon rest mass M#, we refer to the 
ADM mass and the angular momentum of the system, which are given by 



Me 



1 

"27T 



ditjjdS, 

> .... f!_ l AijAV - -K 2 - R k e- 4 



-I 



S^ + l-fA/d^-^A^ k d k f 



d A x, 



d 3 3 



(72) 



(73) 



where dSj = r 2 djrd(cos6)d(p, ^ — —y(d x y + x(d y y , and ip = . In this paper, simulations are performed in axial 
symmetry, and hence, J is conserved. M is approximately conserved since the emission of gravitational waves is 
negligible. Thus, conservation of these quantities is checked during numerical simulations. 
The violation of the Hamiltonian constraint is locally measured by the equation as 



*P 5 



U= Aip - ^R k + 27rp H ^ 5 + \- (An A* - -K 2 



1^6 fei 



tp 5 / ~ 2 



|AV| + If R k k \ + |27rpHV & l + ^-[A ij A» + -K' 



Following [34] , we define and monitor a global quantity as 



H 



m7 / p * Ud3x - 



(74) 



(75) 



Hereafter, this quantity will be referred to as the averaged violation of the Hamiltonian constraint. 



III. NUMERICAL SCHEME FOR SOLVING GRMHD EQUATIONS 



A. GRMHD equations 



As described in Sec. II C, we write the GRMHD equations in the conservative form. In this case, roughly speaking, 
there are two options for numerically handling the transport terms [38]. One is to use the Godunov-type, approximate 
Riemann solver [39,40], and the other is to use the high-resolution central (HRC) scheme [41,20]. We adopt a HRC 
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scheme proposed by Kurganov and Tadmor [42] and very recently used in special relativistic simulations by Lucas- 
Sarrano et al. [43] . Thus our numerical scheme for a solution of GRMHD equations is slightly different from that in 
[10], in which the HLL scheme [44] is basically adopted. 
The basic equations can be schematically written as 

dV , , 91n^ pi = S5 



dt dx l dx i 



where 



\J=(ft*,S i ,So,B i ), (77) 
F' = (p*v j ,S iV j + ae^PtotV - A + e 60 ^tot(^ + P J ) - r% j ,BV - ffV), (78) 

and S denotes the terms associated with the gravitational force. Here, r^f J denotes a magnetic stress defined by 

U j = ^^[B i + Ui (B k u k )}, (79) 

-V = ^r^w . ( 8 °) 

In addition to U, we define a set of variables as 

P = (p„u i ,e,B i ). (81) 

iii and £ are computed at each time step from Eqs. (60) and (62). We use P for the reconstruction of F at cell 
interfaces. In standard method, one often uses a set of primitive variables (p, v % , e, B l ) instead of P for reconstruction 
of F. We have found in the test problems that even using p* and instead of p and v l , it is possible to guarantee 
the similar accuracy and stability. 

To evaluate F, we use a HRC scheme [43]. The fluxes are defined at cell faces. A piece- wise parabolic interpolation 
from the cell centers gives Pr and P^, the primitive variables at the right- and left-hand side of each cell interface, 
as 

QL=Qi+ ^tp^ + m^, (82) 

gs = 0i+ ,_£(r£^_£(5±i^±i (83) 

Here, Q denotes a component of P and A i+1 = Qi + \ — Qi- $ denotes a limiter function defined by 

$(r) = minmod(l, br) (1 < b < 4 for TVD condition), (84) 

where rf = Aj±i/Aj, and 

f 1 if x > 1 

minmod(l,a;) = < x if 1> x>Q . (85) 
[ if x < 

For the simulations presented in Sees. V and VI, we choose 6 = 2 unless otherwise stated. We have found that the 
dissipation is relatively large for b = 1 with which it is difficult to evolve isolated neutron stars for a long time scale 
accurately. On the other hand, for b > 3, the dissipation is so small that instabilities often occur around strong 
discontinuities, and around the region for which P t ot 3> P- 

From Pi and P#, we calculate the maximum wave speed cl and cr, and the fluxes F^ and F^ at the right- and 
left-hand sides of each cell interface. Then, we define c max = max(c,L, cr), and the flux 



F = i 

2 



F L + F R - c max (Uij - U L ) 



(86) 



In adopting the central schemes, the eigen vectors for the Jacobi matrix dF/dXJ are not required in contrast to 
the case with the Godunov-type scheme [38,45]. However, the eigen values for each direction are still necessary to 
evaluate characteristic wave speeds cl and cr. The equation for the seven eigen values A is derived by Anile and 
Pennisi [46]: Three of the seven solutions for A in x l direction are described by 
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A = v\ 



b l ±U l y/ph + b 2 



b^u^ph + b 2 

and rest four are given by the solutions for the following fourth order equation 



( U *)*(A-«T(1-C) + 



( u *) 2 (A-^) 2 (y 



p/i + fe 2 a 
Here, (, the sound velocity c s , and the Alfven velocity v A arc defined, respectively, by 



= (no summation for i). 







2 2 
- v A c s, 




1 


dP 


P dP 




~ h 


[dp 


e p 2 de 


p. 



ph + b 2 ' 



(87) 

(88) 

(89) 
(90) 

(91) 



In the central schemes, we only need the maximum characteristic speed, and thus, only the solutions for Eq. (88), which 
contain the fast mode, are relevant. The solutions for the fourth-order equation are determined either analytically or 
by standard numerical methods. However, for simplicity and for saving computational time, we use the prescription 
proposed by Gammie et al. [20], who have found it convenient to replace the fourth-order equation approximately by 
a second-order one: 



(u 4 ) 2 (A-t/) 2 (l-C)-C( 



r 



(no summation for i). 



(92) 



The solution of Eq. (92) for an arbitrary direction x l is written as 
1 



A* = 



a 2 - V k V k C 



v l a 2 (l-C)- (3\(a 2 -V 2 ) 
± a\/CV(a 2 ~ V 2 ){-f H (a 2 - V 2 () - (1 - QVW} 



(no summation for i), 



(93) 



where V 1 — v l + (i % and V 2 = -fijV l V^ . This is equivalent to that obtained by replacing c s by ^/C m the solution for 
the pure hydrodynamic case [47,48,34]. 



B. Induction equation 

The induction equation may be solved using the same scheme as in solving the hydrodynamic equations described 
above. However, with such a scheme, the violation of the constraint equation (66) is often accumulated with time, 
resulting in a nonrcliablc solution. Thus, we adopt a constraint transport scheme [49]. Namely, we put the com- 
ponents of the magnetic field at the cell-face centers. Here, we specifically consider the axisymmetric case with the 
cylindrical coordinates {x, ip, z) (w is replaced by x). Extension to the nonaxisymmetric case is straightforward, and 
the description below can be used with slight modification. 

In the axisymmetric case with the cylindrical coordinates, the numerical computation is performed in a discretized 
cell for (x, z). Here, we denote the cell center for (x, z) by Then, we put B x at (i + 1/2, j), and B z at + 1/2) 

while components of the gravitational field and fluid variables as well as B v = xB^ are put at the cell center In 
this case, the induction equations for B x and B z are solved in a constraint transport scheme [49] , while that for B v is 
solved in the same method as that for the continuity equation of . 

Computing the flux at cell edges for the induction equation is different from that for the fluid equation. This is 
because numerical fluxes have to be defined so that the constraint equation (66) is satisfied. For example, for the x 
component of the induction equation, the flux in the z direction is written as v z B x — v x B z = F\. On the other hand, 
for the z component of the induction equation, the flux in the x direction is written as v x B z — v z B x = F 2 . Both i*\ 
and F 2 have to be defined at cell edges (i + 1/2, j + 1/2), and for the constraint equation (66) to be satisfied, we have 
to require F\ = —F 2 = F at each cell edge. In addition, an upwind scheme should be adopted for numerical stability: 
For the induction equation of B x , the upwind prescription should be applied for the z component of the flux. On the 
other hand, for the induction equation of B z , the upwind prescription should be applied for the x component of the 
flux. F has to be determined taking into account these requirements. 
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We here adopt a scheme proposed by Del Zanna et al. [41], which satisfies such requirements. In this scheme, the 
flux is written as 

1 r z r x 

p = ^(F LL + F LR + F RL + F RR ) - ™ ax [{B X ) R - (B X ) L ] + ™ ax [(B Z ) R - [B Z ) L ], (94) 

where, e.g., F LR is the flux defined at the left-hand side in the x direction and at right-hand side in the z direction. 
These fluxes are computed by a piece- wise parabolic interpolation. c l max is the characteristic speed for the prescription 
of an upwind flux-construction and calculated at cell edges using the interpolated variables. For simplicity, we set 

c max = max K."fl) and C = max ("L,Wfl)- 

For solving other equations, it is necessary to define the magnetic field at the cell center. Since the x and z 
components of the magnetic field are defined at the cell face centers (i.e., B x at (i + 1/2, j) and B z at + 1/2)), 
this is done by a simple averaging as 

BZ j = l(B? +1/2J +BU/2, j ), (95) 
Bt j = l(Bl j+1/2 + Bl j _ 1/2 ). (96) 

Also, v l at the cell face center is necessary for computing c^ ax and c^ ax . To compute them, we also use a simple 
averaging. For the definition of w 4 fe +1//2 j and v^j +1 ^ 2 , we have also tried the Roe-type averaging in terms of p X J 2 , but 
any significant modification in the results has not been found. 

Before closing this section, we note that our scheme for the induction equation is different from that adopted in 
[10], in which a Toth's method is used [50]. 



IV. INITIAL VALUE PROBLEM 



In the fully general relativistic and dynamical simulations, we have to solve the constraint equations of general 
relativity for preparing the initial condition. One solid method is to give an equilibrium state. For rigidly rotating 
stars of poloidal magnetic fields in axial symmetry, such equilibrium has been already computed [51]. However, for 
differentially rotating stars or nonaxisymmctric cases, the method has not been established. Thus, we here present 
a simple method for preparing an initial condition which is similar to that in [52]. In the following, we assume that 
axisymmetric matter fields /)*, e = hw — P/pw, h, and Uj = hui are a priori given (e.g., those for rotating stars of no 
magnetic field in equilibrium are given). Although we assume the axial symmetry, the same method can be applied 
for the nonaxisymmctric case. 

Initial conditions for magnetic fields have to satisfy Eq. (66). A solution of Eq. (66) is written as 

B k = e kii d i A j , (97) 

where Aj is an arbitrary vector potential and e kl ^ is a Levi-Civita tensor of flat three-space. If we choose 

A X = A Z = 0, and A v ^ 0, (98) 

the magnetic fields are poloidal. Here, we assume to use the cylindrical coordinates (x, ip, z) (w is replaced by x). 
In the axisymmetric case, we can also choose pure toroidal magnetic fields as 

B x = B z = 0, and B v £ 0, (99) 

where B v may be an arbitrary function. In the following, we give a nonzero function either for A v or for B'^. 

Initial conditions also have to satisfy Eqs. (28) and (29). In the following, we assume that jij and K are given 
functions in these equations. Remind that /?h and Ji are written as 

PH = + r 12 (b 2 B2+ Z Ul? ) ' ( 10 °) 

Ji = P . t u t 4>- 6 + ^2 (b 2 u, - BiB j u } , (101) 



where ip denotes the conformal factor (= e^), and w = y/T+ip 4 j i: >UiUj. Thus, if p*, e, h, Uj, A v , B v , K, and 7^ 
are given, the remaining unknown functions are ip and Ajj . This implies that the constraint equations arc solved for 
these variables using the technique developed by York [53] . 
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First, we decompose the tracefree part of the extrinsic curvature as 

An = ^M, = b % w 3 + b w t - ^%b k w k + k^ t , (102) 

where W, is a three vector, and K^ T is a transverse-tracefree tensor which satisfies 

WKj? = Q = Kj?^. (103) 

K^ T would be composed mainly of gravitational waves. Hereafter, we set K^ T = for simplicity. Using Eq. (102), 
Eq. (29) is rewritten to 

AWj + \bjbiW 1 + RjiW 1 - \ip & D K = 87rJ 4 V> 6 . (104) 

This equation can be solved for an initial trial function of ip. Then, Aij is computed from Eq. (102). Substituting 
Aij, the Hamiltonian constraint (28) is solved in the next step. Then we solve the momentum constraint again, and 
repeat these procedures until a sufficient convergence is achieved. 

V. SPECIAL RELATIVISTIC TESTS 

In this section, we present numerical results for a number of special relativistic tests. In the tests, we adopt the 
T-law equations of state as 

P=(T- l)pe, (105) 
with T = 4/3 or 5/3. Simulations are always performed using the uniform grid in all the axis directions. 

A. One dimensional tests 

Any numerical implementation of the MHD equations has to be checked if it can produce the basic waves such 
as shock and rarefaction waves accurately. Komissarov [39] has proposed a suite of one-dimensional test problems 
in special relativity: Propagation of fast and slow shocks, fast and slow rarefaction waves, Alfven waves, compound 
waves, shock tube tests, and collision of two flows. We have performed all the tests except for the compound wave 
following [10]. Our implementation can integrate each of remaining eight tests although in some cases we have to 
reduce the Courant number significantly to avoid numerical instabilities as reported by Gammie [20]. On the other 
hand, we adopt the same limiter, 6 — 2, for all the simulations. Numerical results are shown in Figs. 1-5. Grid size 
N and spacing Ax we adopt for each of test simulations are approximately the same as those by Komissarov, and 
described in the figure captions. 

Figure 1 shows the results for fast and slow shocks. In these problems, the system is stationary with respect to the 
frame comoving with the shock front. The velocity of the shocks is 0.2 and 0.5 for the fast and slow shocks, respectively. 
As the previous works illustrate [39,20,21,10], the fast shock can be computed accurately with a relatively large grid 
spacing. On the other hand, in the numerical solution of the slow shock, a spurious modulation is found for p in the 
region of 1 <, x <, 1.3 as in the previous works [39,20,21,10]. This is always generated soon after the onset of the 
simulation irrespective of grid resolutions. Thus, it is impossible to avoid such small error in our implementation. 
Although the modulation is always present, its wavelength and amplitude gradually decrease with improving the 
grid resolution. We computed an LI norm defined for the difference between the numerical and exact solutions, and 
found that it decreases as the grid spacing is smaller. In this case, the convergence is achieved at first order since 
discontinuities are present, around which the transport terms of hydrodynamic equations are computed with the 
first-order accuracy. 

Figure 2 shows the results for switch-off and switch-on rarefaction waves. Although we have not compared the 
results precisely with those by other authors [39,20,21,10], the accuracy of our results is similar to that reported by 
others. For the switch-off waves, a spurious bump is found at x ~ —0.6 as in the previous works [39,20,21,10]. As 
in the slow shock problem, this bump is generated at t = irrespective of grid resolutions, and with improving the 
grid resolution, the magnitude of the LI norm decreases at first order. On the other hand, the numerical solution for 
the switch-on waves, spurious bumps are not present, and with Ax = 0.005, a good convergent result appears to be 
obtained with our implementation. 
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FIG. 1. (a) Propagation of a fast shock. The snapshot at t = 2.5 is shown. Numerical simulation is performed with N = 100 
and Ax = 0.02. (b) Propagation of a slow shock. The snapshot at t = 2 is shown. The numerical simulation is performed with 
N = 400 and Ax — 0.005. Only one fourth of data points are plotted except for [0.92,1.08] in which all the data points are 
plotted. For both cases, the initial discontinuities were located at x = 0, and the shock fronts move with a constant velocity \i 
where \i = 0.2 and 0.5 for the fast and slow shocks, respectively. 




(a) (b) 

FIG. 2. (a) Propagation of a fast switch-off rarefaction waves. The snapshot at t = 1.0 is shown, (b) Propagation of a slow 
switch-on rarefaction waves. The snapshot at t = 2 is shown. For both cases, the initial discontinuities were located at x = 0, 
and the numerical simulations are performed with N = 400 and Ax = 0.005. Only half of all the data are plotted for both 
figures. 
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FIG. 3. (a) Propagation of a strong continuous Alfven waves. The snapshot at t = 2 is shown. The initial configurations 
are the analytic solution at t = 2 are shown by the dashed curves. The numerical simulation is performed with N = 500 and 
Aa; = 0.005. Only 1/4 of all the data are plotted, (b) LI norm of the error for p (circles) and P (triangles) as a function of 
Ax. These variables should be constant in this problem, and the deviation from the stationary values is due to a numerical 
error. The dotted line denotes the slope expected for second-order convergence. 




(a) (b) 

FIG. 4. (a) Shock tube problem with the initial discontinuities at x — normal to the magnetic field. The numerical 
simulation is performed with N — 500 and Aa; = 0.005. We note that a large bump at x ~ 0.9 for u x is due to a numerical 
error associated with the limiter (6 = 2) of a weak dissipation. The height of this bump is decreased if we use a more dissipative 
limiter (e.g., the minmod limiter with b = 1). (b) Shock tube problem with the initial discontinuities at x — parallel to the 
magnetic field. The numerical simulation is performed with N = 500 and Aa; = 0.005. For both cases, the snapshot at t = 1.0 
is shown. 
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FIG. 5. Collision of two flows with opposite directions of the tangential component of the magnetic field. The snapshot at 
t = 1.22 is shown. The numerical simulation is performed with N = 400 and Ax — 0.005. 

Figure 3(a) shows the results for an Alfvcn wave test, demonstrating that the Alfven wave can be computed 
accurately with our implementation as in [39,10]. In this problem, the density and pressure should be unchanged. In 
our results, this is achieved within ~ 1% error for Ax — 0.0025. Since no discontinuities are present in this problem, 
the convergence of the numerical solution to the exact one should be achieved approximately at second order [10] . To 
check if this is the case, we compute an LI norm defined by the difference between the numerical and exact solutions 
for p and P. The results are shown in Fig. 3(b), which illustrates that the convergence is achieved approximately at 
second order (slightly better than second order). 

In Fig. 4, numerical results for shock-tube problems are presented. For the problem of Fig. 4(a), shocks are very 
strong since the ratio of the pressure in the left- and right-hand sides at t = is 10 3 . However, since the magnetic 
field lines are normal to the discontinuities, the effects of the magnetic field for the formation and propagation of 
shocks are absent. In this case, a large spurious overshooting is found around the shock for u x . This is partly due 
to our limiter (b = 2) which is not very dissipative. If we use the minmod limiter (b = 1), height of the overshooting 
decreases although the shocks are less sharply computed. For the problem of Fig. 4(b), shocks are not as strong as 
those in 4(a). However, the magnetic fields affect the formation and propagation of shocks since they are parallel 
to the shocks. The results shown in Fig. 4(b) are very similar to those in [39,20,21,10], and hence, are likely to be 
as accurate as them. This indicates that our implementation can compute magnetized shocks as accurately as the 
previous ones. 

In Fig. 5, numerical results for collision of two magnetized flows are presented. It shows that four separate 
discontinuities generated at t = are computed accurately. As found in previous papers [20,10], a small dip spuriously 
appears in p around x = 0. As in the case of the slow shock and switch-off rarefaction wave, this is spuriously generated 
at t = irrespective of grid resolution, and with improving the resolution, the magnitude of the error is decreased at 
first order. 



For multidimensional tests, following Del Zanna et al. [41], we performed simulations for (i) a cylindrical blast 
explosion, (ii) a rotating cylinder in two-dimensional Cartesian coordinates with a uniform magnetized medium, 
and (iii) propagation of a jet in cylindrical coordinates in a magnetized background. The parameters for the initial 
conditions adopted here are the same as those in [41]. On the other hand, we varied the grid spacing to see the 
convergence in contrast to the previous works. 

In the test (i), the Cartesian grid of (x, y) is adopted with the range [—0.6, 0.6] for both directions. The grid spacing 
chosen is 0.004, 0.005, and 0.006. The initial condition is 



B. Multi dimensional tests 




(106) 



15 




FIG. 6. Snapshot of cylindrical blast explosion (multidimensional test (i)) at t — 0.4. The contour curves for the density, 
pressure, and magnetic pressure P mag , as well as the magnetic field lines are shown. The contour curves for each quantity 
(denoted by Q) are drawn for Q = Qmax x io~°' 4xi (i — 1-8) (solid curves) and Q = Q m ax x 10 -0 ' 01 (thick solid curves). Here 
Qmax denotes the maximum value. The results with Aa; = 0.004 are presented. 
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FIG. 8. The same as Fig. 6 but for an explosion of a rotating cylinder (multidimensional test (ii)) at t = 0.4. The grid 
spacing for the corresponding simulation is 0.0025. 
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FIG. 9. Configuration of various quantities in an explosion of a rotating cylinder at t = 0.4 with different grid resolutions 
(A denotes the grid spacing). 
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FIG. 10. The density contour curves and velocity vectors (left panel) and the magnetic field lines (right panel) at t = 15 
and t — 30 for an axisymmetric jet. The grid spacing for the corresponding simulation is 0.06. The density contour curves are 
drawn by the same method as in Fig. 6. 



with u % = and T = 4/3. Because of the large internal energy in the central region, the outward explosion occurs. In 
this problem, the shocks generated at t = are very strong, and hence, the minmod limiter with b = 1 is adopted to 
avoid numerical instability. With the limiter of b — 2, the computation crashes because of the appearance of negative 
internal energy (or h < 1) irrespective of the grid resolution. We have found that the computation first crashes along 
the line of x = y and x = —y for which the accuracy is likely to be worst. 

In Fig. 6, we display the snapshot of the numerical results at t = 0.4. In Fig. 7, configurations of the density, 
pressure, magnetic pressure, and Lorentz factor along x and y axes are shown for three levels of grid resolution. The 
expansion velocity of the blast wave is largest along the x axis because of the confinement by the magnetic pressure. 
The maximum value of the Lorentz factor is about 4 at t = 0.4 with the best resolved case. Along the y axis, the 
magnetic field lines arc squeezed yielding the highest magnetic pressure. These features agree with those found in 
[39,41]. As mentioned in [41], the total energy is completely conserved since we solve the MHD equations in the 
conservative form and do not add any dissipative terms in contrast with the treatment in [39] . 

One point to be mentioned is that convergence around the density peak along the x axis is not achieved well 
within the adopted resolution although for other region, convergence is achieved well. The likely reason is that the 
discontinuities around the peak is very thin for which it is very difficult to resolve with the chosen grid resolutions. 
Thus, it is difficult to accurately derive the maximum values of the density, pressure, and Lorentz factor which are 
underestimated in this test problem. 

In the test (ii), the Cartesian grid of (x,y) is also adopted with the range [—0.6,0.6] for both directions. The grid 
spacing chosen is 0.0025, 0.003, and 0.004. The initial condition is 



(o pb x m-( ( 10 ' 1 ' 1 '°) for v2±Z<o.i (W7) 

(P ^^ '° j ~ I (1,1,1,0) for V?+?>0.1, 1 ' 

with v 1 = [— uiy,ujx] for \/ x 2 + y 2 < 0.1 where lu = 0.995 and thus the Lorentz factor at the surface of the rotating 
cylinder is initially about 10. T is chosen to be 5/3 following [41]. 

In Fig. 8, we display the snapshot of the numerical results at t — 0.4. In Fig. 9, configurations of the density, 
pressure, magnetic pressure, and Lorentz factor along x and y axes are shown for three levels of grid resolution. In 
this problem, the magnetic field lines keep winding-up, and at t = 0.4, the central field lines are rotated by an angle 
of ~ 90 degrees. Because of magnetic braking, the rotational speed is decreased monotonically, and at t = 0.4 the 
maximum Lorentz factor is decreased to ~ 1.7. Due to the outward explosion induced by the rotation, the density in 
the central region becomes an uniformly low value of ~ 0.44 while an ellipsoidal density peak is formed around the 
central region. As in the test (i), it is difficult to obtain a convergent value for the peak density with the chosen grid 
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resolutions. The likely reason is that the thickness of the density peak is so small that the grid resolutions are not 
sufficient. However, for the other region, convergent results are obtained. 

In the test (iii), the cylindrical grid of (x,z) is adopted with the range [0,8] and [0,20], respectively. The grid 
spacing is 0.06, 0.08, and 0.1. The initial condition is 

- { arw — r 1 - 1 m 

with v x = 0, B x = 0, and T = 5/3. The region with < x < 1 and < z < 1 is defined to be a jet-inlet zone, and the 
stationary condition is artificially imposed. In the simulation, the regularity condition is imposed along the symmetric 
axis x = 0. For the boundary conditions at z = 0, extrapolation is assumed following [41]. In this test, we adopt the 
minmod limiter with b — 1 since with 6 = 2, the computation soon crashes irrespective of the grid resolutions. 

In Fig. 10, we show the snapshot of the density contour curves and magnetic field lines at t = 15 and 30 with 
Ax = 0.06. The contour curves and field lines are similar to those in [41]. The maximum Lorentz factor at t = 
is <~ 7.09. At the head of the jet, the density becomes maximum and shocks are formed, inducing back flows at 
the shocks. These flows make a cocoon which is to expand in the direction of the cylindrical radius, squeezing the 
magnetic field lines. A part of the matter is back-scattered toward the z = plane dragging the magnetic field lines 
together. As a result, the magnetic field lines are highly deformed. The deformation is computed more accurately 
with finer grid resolutions. 

However, we found that precise computation for the deformation of the magnetic field lines increases the risk for 
crash of the computation. For Ace = 0.1 and 0.08, computations can be continued until the shock front of the jets 
reaches the outer boundary. However, for Ax = 0.06, the computation crashes at t ~ 35 in spite of the fact that the 
motion of the jet head is still stably computed. If we adopt a better resolution with Ax ^ 0.05, the computation 
crashes before t reaches 30. The instabilities always occur near the boundary region of the jet-inlet zone around which 
the magnetic field configuration is deformed to be highly complicated. This seems to be due to the fact that we 
impose the stationary condition inside the jet-inlet zone. This artificial handling makes the field configuration near 
the boundary of the jet-inlet zone nonsmooth (i.e., the derivative of the magnetic field variables can be artificially 
larger for better grid resolutions). Here, we note that this problem happens only in the presence of magnetic fields. 
Thus, for continuing the computation for a longer time, probably, it is necessary to include a resistivity for inducing 
reconnections of magnetic fields near the jet-inlet zone for stabilization. The other method is to change the stationary 
condition wc adopt here to other appropriate boundary conditions near the jet- inlet zone [54]. 

VI. GENERAL RELATIVISTIC TESTS 

A. Relativistic Bondi accretion 

As the first test for general relativistic implementation, we perform a simulation for spherical accretion onto the 
fixed background of a Schwarzschild black hole. The relativistic Bondi solution is known to describe a stationary flow, 
and thus, by comparing the numerical solution with the analytical one, it is possible to check the suitability of the 
numerical implementation for general relativistic hydrodynamics problems [55,56]. Furthermore, it has been shown 
that the relativistic Bondi solution is unchanged even in the presence of a divergence-free pure radial magnetic field 
[21]. Thus, it can be also used for checking the GRMHD implementations. The advantage of this test is that the 
exact solution can be obtained very easily while it involves strong gravitational fields, relativistic flows, and strong 
magnetic fields all together. 

Following previous authors [20,10], we write the metric in Kerr-Schild coordinates in which all the variables are 
well behaved at the event horizon (r = 2M; where r and M are the radial coordinate and the mass of the black 
hole). Nevertheless, the hydrostatic equations for the stationary flow are the same forms as those in the Schwarzschild 
coordinates, and thus, the stationary solution is determined from an algebraic equation which can be easily solved by 
standard numerical methods [57]. 

For this test, we adopt the same solution used in [55,20,21,10]. Namely, the sonic radius is set at r = 8M, the 
accretion rate M = 4irpr 2 u r is set to be —1, and the adiabatic index for the equation of state is 4/3. The simulation 
is performed in an axisymmctric implementation with the cylindrical coordinates (x, z). The computational domain 
is set to be [0, 18M] for x and z, and the radius of r = 1.9M is chosen as the excision radius. The uniform grid is 
adopted. At the excision radius and outer boundaries, we impose the condition that the system is stationary. The 
(semi) analytic solution for the stationary Bondi flow is put as the initial condition, and we evolve for 100M following 
previous authors [20,21,10]. 
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The simulations are performed changing the grid spacing Ax. Irrespective of the grid resolution, the system relaxes 
to a stationary state long before 100M. When evolved with a finite-difference implementation, discretization errors 
will cause small deviations in the flow from the exact stationary configuration. These deviations should converge 
to zero at second order with improving the grid resolution. To diagnose the behavior of our numerical solution, we 
measure an LI norm for — p° xact where p® xact denotes the exact stationary value of p*. Specifically, the LI norm 
is here defined by 

/ \p,^ P r ct \d 3 x/ f P T act d 3 x. (109) 

Jr>2M 1 Jr>1M 

For the convergence test, the grid spacing is changed from 0.06M to QAM . 

The radial magnetic field strength is also changed for a wide range. Following [20,10], we denote the magnetic field 
strength by 



P 



r=2M 



(110) 



We note that the ratio of the magnetic pressure to the gas pressure b 2 /2P is w 3.85/3 at r = 2M for the solution 
chosen in this test problem. 
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FIG. 11. LI norm for the error of density p* as a function of the grid spacing in units of M for the magnetized Bondi flow. 
The numerical numbers attached for each curve denote the values of (3 (see Eq. (110) for its definition). 



In Fig. 11, we show the LI norm as a function of the grid spacing for < < 63. Irrespective of the magnetic field 
strength, the numerical solution converges approximately at second order to the exact solution for Ax — > 0. The LI 
norm is larger for the stronger magnetic fields, implying that the relaxed state deviates more from the true stationary 
solution for the larger value of 0. Specifically, the velocity field configuration deviates significantly from the exact 
solution with increasing the value of 0, although the deviation for the density configuration is not very outstanding. 

In this test simulation, we have found several interesting behaviors of the numerical solutions. First, for a given 
value of the grid spacing with Ax <; 0.1M, there is the maximum allowed value of above which the computation 
crashes. The maximum value is larger for better grid resolution; e.g. for Ax — 0.1M, 0.2M, and 0.3M, the maximum 
allowed values of are w 45, 25, and 10, respectively. For Ax ^ 0.1M, on the other hand, the maximum allowed value 
is <~ 70 irrespective of the grid resolution. The limitation is due to the well-known weak point in the conservative 
scheme that the small error in the magnetic energy density in the magnetically dominated region with ^> 1 leads to 
fractionally large errors in other components of the total energy density, by which the computation crashes (typically, 
the internal energy density becomes negative). For the poorer grid resolutions, the numerical error is larger, and hence, 
the computation crashes for the smaller value of the magnetic field. Second, the maximum allowed value of found 
here (~ 70) is by about one order of magnitude smaller than that found in [20]. This is probably due to the difference 
of the coordinate system adopted; we use the cylindrical coordinates while the authors in [20] use the spherical polar 
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coordinates which obviously have advantage for handling the spherically symmetric problem. However, we note that 
even in the cylindrical coordinates, it is possible to handle the flow with a very high value of ~ 60 if a sufficient 
grid resolution is guaranteed. In [10], the authors suggest that in the cylindrical coordinates, the maximum allowed 
value of is at most ~ 10. We have not found such severe limitation in our numerical experiment. Their failure 
for simulating the flow with high values of is probably due to the fact that they use an excision boundary which 
may be applicable for more general problems (e.g., for simulation of dynamical spacetimes). Even in the cylindrical 
coordinates, a high value of will be achieved if a stationary inner boundary condition is imposed. 



B. Longterm evolution for system of a rotating star and a disk with no magnetic field 

Next, we illustrate that with our implementation (for axisymmetric systems), self-gravitating objects can be simu- 
lated accurately. In a previous paper [61], we have already illustrated that our implementation with a HRC scheme 
can simulate rapidly rotating compact neutron stars for more than 20 rotational periods accurately. Thus, we here 
choose a more complicated system; an equilibrium system composed of a rapidly rotating neutron star and a massive 
disk. By this test, it is possible to check that our implementation is applicable to a longterm evolution not only for 
an isolated rotating star but also for a self-gravitating disk rotating around a compact object. 

The equilibrium configuration is determined by solving equations for the gravitational field and hydrostatic equations 
self consistently. For simplicity, we here adopt a conformally flat formalism for the spatial metric [58] . As shown in 
[59], a good approximate solution for axisymmetric rotating stars can be obtained even in this approximation. Thus, 
the initial condition presented here can be regarded as a slightly perturbed equilibrium state. At the start of the 
simulations, we further add a slight perturbation by reducing the pressure by 0.1% to investigate if a quasiradial 
oscillation is followed stably and accurately. The magnitude of the perturbation in association with the conformally 
flat approximation is much smaller than this pressure perturbation. 

The Euler equation for axisymmetric stars in equilibrium can be integrated to give the first integral, which is written 

as 

ln\ + [ v^u^dQ, = C, (111) 
u l J 

or 

+ J hu v dn = C, (112) 

where C and C are integral constants. Equation (111) is a well-known form [60]. However here, we adopt Eq. (112), 
and set that the specific angular momentum hu v is constant (= j) for the disk and O =const for the central star. 

A hybrid, parametric equation of state is used in this simulation following previous papers [15,26,8] . In this equation 
of state, one assumes that the pressure consists of the sum of polytropic and thermal parts as 

P = P P + P th . (113) 

The polytropic part, which denotes the cold part of the equations of state, is given by 

P? = \ K K lpT r; P i Pmc ' (114) 

[ K 2 p l2 , p > /9nuc, 

where K\ and K 2 are polytropic constants. p nuc denotes the nuclear density and is set to be 2 x 10 14 g/cm 3 . In this 
paper, we choose I\ = 4/3 and T 2 = 2.5. Since Pp should be continuous, the relation, K 2 = Kip^^ 2 > 1S required. 
Here, the value of K\ is chosen to be 2.5534 x 10 14 in the cgs unit. With this value, the maximum ADM (baryon 
rest) mass for the cold and spherical neutron star becomes about 1.84M© (2.05M Q ) which is a close to that derived 
in realistic equations of state [62] . 

Since the specific internal energy should be also continuous at p = p nuc , the polytropic specific internal energy ep 
is defined as 



£p = < 



Ki r i 



Vr 2 _, ( 115 ) 
I r 2 - i + (iwxrw) ' p^p^- 
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The thermal part of the pressure P t h plays an important role in the case that shocks are generated. P t h is related to 
the thermal energy density e t h = £ — £p as 

Ph = (Tth - l) P s th . (116) 

For simplicity, the value of r t h, which determines the strength of shocks, is chosen to be equal to Ti. For computing 
initial equilibria, we set e = sp and P = Pp. 
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(a) (b) 
FIG. 12. (a) Evolution of the central density of a neutron star and mass and angular momentum of a disk around the neutron 
star, (b) Evolution of the ADM mass, angular momentum, and averaged violation of the Hamiltonian constraint. The ADM 
mass and angular momentum are shown in units of their initial values. The solid, dashed, and long-dashed curves show the 
results with (241,241), (193, 193), and (161,161) grid sizes, respectively. 




For the simulation, we choose a sufficiently deformed star with the axial ratio of the minor axis to major axis 
~ 0.6. The ADM mass is 1.888M , total baryon rest mass 2.074Af Q , the central density 1.3 x 10 15 g/cm 3 , the 
circumferential radius at equator 16.2 km, the rotational period P c = 0.841 ms, and J/M 2 = 0.545. Thus, the 
neutron star considered is massive and rapidly rotating. The baryon rest mass of the disk is much smaller than the 
central star as 4.9 x 10~ 5 M Q with the maximum density s=s 2 x 10 10 g/cm 3 . Since it is of low density, the disk is 
composed of r = 4/3 polytropic equation of state. Orbital radius of inner edges of the disk is ~ 20 km, and thus, the 
uniform specific angular momentum is small as j w 3.45M which is very close to the value for a particle orbiting an 
innermost stable circular orbit. The rotational periods of the disk at the inner and outer edges in the equatorial plane 
are 1.03 ms(= 1.2P C ) and 2.58 ms(= 3.1P C ), respectively. The simulations are performed in axial symmetry with 
(241,241), (193, 193), and (161,161) grid sizes for which the grid spacing is 0.165, 0.202, and 0.248 km, respectively. 
The reflection symmetry with respect to the equatorial plane is assumed. The outer boundaries along each axis are 
located at 39.6 km. 

An atmosphere of small density p = 2 x 10 4 g/cm 3 is added uniformly outside the neutron star and disk at 
t = 0, since the vacuum is not allowed in grid-based hydrodynamics implementations. We note that the density of 
atmosphere can be chosen to be much smaller than the nuclear density p nuc - This is the advantage of HRC schemes 
in which such low density can be handled in contrast with high-resolution shock-capturing schemes [61]. Since the 
atmosphere is added as well as a small pressure perturbation is imposed, the Hamiltonian and momentum constraints 
are enforced at t = using the method described in Sec. IV. 

In Fig. 12(a), we show the evolution of the central density of the neutron star, and mass and angular momentum 
of the disk which are defined by 



M, 



J X 

=-j 



x. 



(117) 



(118) 



where x- m denotes the initial coordinate radius of the inner edges of the disk. The figure shows that our implementation 
keeps the equilibrium system to be in equilibrium for more than 20P C . With the grid of size (241,241), increase of 
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the density, which is perhaps associated with the outward transport of the angular momentum, is at most ~ 1% at 
t = 20P C . The change in the baryon rest mass and angular momentum of the disk, which is caused spuriously by the 
mass transfer from the central star and mass accretion to the central star due to a numerical error, is smaller than 
<~ 0.1%. Also, the numerical results converge at better than second order with improving the grid resolution. 

In Fig. 12(b), we also show the evolution of the ADM mass, angular momentum, and averaged violation of the 
Hamiltonian constraint. It is found that the ADM mass is conserved within <~ 1% error for t < 20P C with (241, 241) 
grid resolution. An outstanding feature is that the angular momentum is conserved with much better accuracy than 
the ADM mass. This is a feature when a HRC scheme is adopted [61]. The averaged violation of the Hamiltonian 
constraint also remains to be a small magnitude for t < 20P C and converges at better than second order. All these 
results illustrate that our implementation can compute self-gravitating equilibrium systems accurately. 



Next, we add magnetic fields confined only in the disk around the neutron star. For this test, we use the same 
system of a neutron star and a disk which is described in Sec. VI B. Similar test in a fixed background spacetime 
of a black hole has been performed in [20,21]. Here, we perform the test in full general relativity replacing the black 
hole by a neutron star. The purpose of this subsection is to illustrate that our implementation can follow the growth 
of magnetic fields by winding-up due to differential rotation of the disk. Subsequent papers will focus on detailed 
scientific aspect of this issue [63]. 

Following [20,21], the ip component of the vector potential A v is chosen as 



where A is a constant which determines the magnetic field strength. Then the magnetic fields are given by B z = 
x~ 1 d x A v and B x = —x~ 1 d z A ip . This choice of A v produces poloidal field loops that coincide with isodensity contours. 
Here, po is chosen as 0.3/? max: disk where /c m ax:disk is the maximum density inside the disk. In the following, all the 
simulations are performed in axial symmetry with (301, 301) grid size and with the grid spacing of 0.165 km. The 
reflection symmetry with respect to the equatorial plane is assumed. We note that the boundary condition for the 
magnetic field is B x — B v — and d z B z = at the equatorial plane (in contrast to those for velocity fields v l and Ui 
for which, e.g., d z v x = d z v y = and v z — at the equatorial plane). Outer boundary conditions are not necessary for 
the magnetic field in the present simulations since the location of the outer boundary is far enough from the center 
that the magnetic field lines do not reach the outer boundaries. The Hamiltonian and momentum constraints are 
enforced at t = using the method described in Sec. IV. Since the magnetic field strength we choose is very weak 
initially, the obtained initial condition is approximately the same as that of no magnetic fields. 

Simulations are performed for various values of A which is chosen so that the magnetic pressure is initially much 
smaller than the gas pressure. In the following we specify the model in terms of the initial ratio of the energy of 
magnetic fields to the internal energy of the disk (hereafter Rb) instead of A. Here, the energy of magnetic fields and 
the internal energy of the disk is simply defined by 



and thus, Rb = CAn ag /[/disk at i = 0. We note that the precise definition of the magnetic energy is unknown in general 
relativity, but the present definition is likely to give a guideline for the magnitude within an error of a factor of ~ 2-3. 

In Fig. 13, we show the evolution of U mag for three values of Rb- Here, the magnetic energy is plotted in units of 
the initial value of t/disk (hereafter i7disko) which is about 1.8 x 10 _4 Mdi s k- It is found that f7 mag grows monotonically 

1 /2 

until the growth is saturated irrespective of the value of Rb- The growth rate is in proportional to R B ' in the 
early phase before the saturation is reached. This indicates that differential rotation winds up the magnetic field 
lines for amplifying the field strength [12,13]. After the saturation occurs, i7 mag /[/disko relaxes to <~ 0.02-0.2. These 
values indicate that the magnetic (3 parameter often referred in [21] is of order ~ 10. These relaxed values are in 
good agreement with previous results obtained in the simulation with a fixed background [21]. The magnetic energy 
reached after the saturation depends on Rb, indicating that not only the winding-up of the field lines but also other 
mechanisms (which may be MRI or other instabilities associated with the magnetic fields) are likely to determine the 
final value. 



C. Winding-up of magnetic field lines in a disk around a neutron star 
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FIG. 13. Evolution of [/ mag in unit of initial internal energy of C/disko for Rb = 1.3 x 10 -5 (long-dashed curve), 5 x 10 -5 
(solid curve), and 2 x 10 -4 (dashed curve). In the smaller panel, we display evolution of C/ mag /t/disko for Rb = 5 x 10 -5 with 
three levels of grid resolutions with sizes (301, 301) (solid curve), (241, 241) (dashed curve), and (201, 201) (long-dashed curve). 
Three curves approximately agree. 

To check that the growth of the magnetic fields occurs irrespective of grid resolution, we performed additional 
simulations for Rb — 5 x 10~ 5 with grid sizes of (241, 241) and (201, 201) without changing the location of the outer 
boundaries. In the small panel of Fig. 13, evolution of the magnetic energy for these cases as well as for (301, 301) grid 
size is displayed. It is shown that the growth rate depends very weakly on the grid resolution. This confirms that our 
simulation can follow the winding-up of the magnetic field lines well. On the other hand, it should be mentioned that 
the fastest growing mode of the MRI cannot be resolved in the present computational setting since the characteristic 
wavelength for this mode ~ 2itva/Q, where va denotes the characteristic Alfven speed, is approximately as small as 
the grid size (for Rb = 5 x 10~ 5 , 2ttva/Q ~ Ax) in the current setting. To resolve the fastest growing mode, the grid 
size should be at least one tenth of the present one. Performing such a simulation of high resolution is beyond scope 
of this paper and an issue for the next step. 

In Fig. 14, snapshots of the density profile of disks are displayed for which Rb = 2 x 10~ 4 . It shows that with 
the growth of magnetic fields due to winding-up of the field lines, a wind is induced to blow the matter in the outer 
part of the disk off. Also, the matter in the inner part of the disk gradually falls into the neutron star because of 
the angular momentum transport by the magnetic fields from the inner to the outer parts (see Fig. 15). After the 
nonlinear development of the turbulence, the disk settles down to a quasi stationary state. As explained in [64,21], 
this is probably due to the imposition of axial symmetry which precludes the development of the azimuthal unstable 
modes. Also, in the present numerical simulation, MRI which could induce turbulence is not well resolved. This may 
be also a reason. 

In Fig. 15, we show the evolution of mass and angular momentum of the disk. It shows that after the saturation 
of the nonlinear growth of the magnetic fields, these quantities decrease. Decrease rates of the mass and angular 
momentum take maximum values soon after the growth of the magnetic pressure is saturated (e.g., at t ss 5P C for 
Rb = 2 x 10~ 4 ; cf. the dashed curves). Then, the mass and angular momentum relax to approximately constants 
(cf. the dashed curves). This indicates that the disk settles down to a quasistationary state. An interesting feature 
is that Jdisk is approximately proportional to Md; s k throughout the evolution. This is reasonable because the specific 
angular momentum j is constant in the disk at t = 0, and approximately so is the matter fallen to the neutron star 
as long as the magnetic pressure is much smaller than the gas pressure. However, in the case of Rb = 2 x 10~ 4 , at 
t ~ 20P C for which growth of the magnetic field has already saturated enough, Jdisk /-^disk slightly deviates from the 
initial value. This indicates that angular momentum is transported by the effect of magnetic fields. 

We also performed a simulation for a toroidal magnetic field B v . For B v , we gave 



where po = 0.3/9 max: disk and zq is a constant much smaller than the scale hight of the disk. We note that B v has to be 
zero in the equatorial plane because we impose a reflection symmetry for the matter field with respect to this plane. 
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FIG. 14. Snapshots of the density profile of a neutron and a disk in x — z plane at t = 0, 9.78, and 19.40 ms for Rb = 2 x 10 . 
For p > 10 10 g/cm 3 , the density is denoted by the same color (red), and for 10 10 g/cm 3 > p > 10 7 g/cm 3 , the color is changed 
(from red to green). 
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FIG. 15. Evolution of Mdisk and Jdisk for Rb = (dotted curves), 1.3 x 10 -5 (long-dashed curves), 5 x 10" 5 (solid curves , 
and 2 x 10~ 4 (dashed curves). 

In this simulation, magnetic energy decreases monotonically due to a small expansion of the disk induced by the 
magnetic pressure. In this case, no instability sets in. This is a natural consequence since the field lines are parallel 
to the rotational motion, and hence, they are not wound by the differential rotation. Obviously, the assumption of 
the axial symmetry prohibits deformation of the magnetic field lines and plays a crucial role for stabilization. If a 
nonaxisymmetric simulation is performed, MRI could set in [64,65,23]. 

VII. SUMMARY AND DISCUSSION 

In this paper, we describe our new implementation for ideal GRMHD simulations. In this implementation, Einstein's 
evolution equations are evolved by a latest version of BSSN formalism, the MHD equations by a HRC scheme, and 
the induction equation by a constraint transport method. We performed a number of simulations for standard test 
problems in relativistic MHD including special relativistic magnetized shocks, general relativistic magnetized Bondi 
flow in the stationary spacetime, and fully general relativistic simulation for a self-gravitating system composed of a 
neutron star and a disk. Our implementation yields accurate and convergent results for all these test problems. In 
addition, we performed simulations for a magnetized accretion disk around a neutron star in full general relativity. It 



25 



is shown that magnetic fields in the differentially rotating disk are wound, and as a result, the magnetic field strength 
increases monotonically until a saturation is achieved. This illustrates that our implementation can be applied for 
investigation of growth of magnetic fields in self-gravitating systems. 

In the future, we will perform a wide variety of simulations including magnetized stellar core collapse, MRI for 
self- gravitating neutron stars and disks, magnetic braking of differentially rotating neutron stars, and merger of 
binary magnetized neutron stars. Currently, we consider that the primary target is stellar core collapse of a strongly 
magnetized star to a black hole and a neutron star which could be a central engine of gamma-ray bursts. Recently, 
simulations aiming at clarifying these high energy phenomena have been performed [66,67]. In such simulations, 
however, one neglects self-gravity and also assumes the configuration of the disks around the central compact object 
and magnetic fields without physical reasons. On the other hand, Newtonian MHD simulations including self-gravity 
consistently have recently performed in [68] . However, stellar core collapse to a black hole and gamma-ray bursts are 
relativistic phenomena. For a self consistent study, it is obviously necessary to perform a general relativistic simulation 
from the onset of stellar core collapse throughout formation of a neutron star or a black hole with surrounding disks. 
Subsequent phenomena such as ejection of jets and onset of MRI of disks should be investigated using the output of 
the collapse simulation. In previous papers [26,9], we performed fully general relativistic simulations of stellar core 
collapse to formation of a neutron star and a black hole in the absence of magnetic fields. As an extension of the 
previous work, simulation for stellar core collapse with a strongly magnetized massive star should be a natural next 
target. 

It is also important and interesting to clarify how MRI sets in and how long the time scale for the angular momentum 
transport after the onset of the MRI is in differentially rotating neutron stars. Recent numerical simulations for 
merger of binary neutron stars in full general relativity [6,7] have clarified that if total mass of the system is smaller 
than a critical value, the outcome after the merger will be a hypermassive neutron star for which the self-gravity is 
supported by strong centrifugal force generated by rapid and differential rotation. Furthermore, the latest simulations 
have clarified that the hypermassive neutron star is likely to have an ellipsoidal shape with a large ellipticity [7], 
implying that it can be a strong emitter of high-frequency gravitational waves which may be detected by advanced 
laser interferometric gravitational wave detectors [69] . In our estimation of amplitude of gravitational waves [69] , we 
assume that there is no magnetic field in the neutron stars. However, the neutron stars in nature are magnetized, and 
hence, the hypermassive neutron stars should also be. If the differential rotation of the hypermassive neutron stars 
amplifies the seed magnetic field via winding-up of magnetic fields or MRI very rapidly, the angular momentum may 
be redistributed and hence the structure of the hypermassive neutron stars may be significantly changed. In [7], we 
evaluate the emission time scale of gravitational waves for the hypermassive neutron stars is typically ~ 50-100 ms 
for the mass M ~ 2.4-2.7M Q assuming the absence of the magnetic effects. Here, the time scale of ~ 50-100 ms is 
an approximate dissipation time scale of angular momentum via gravitational radiation, and hence in this case, after 
~ 50-100 ms, the hypermassive neutron stars collapse to a black hole because the centrifugal force is weaken. Thus, it 
is interesting to ask if the dissipation and/or transport time scale of angular momentum by magnetic fields is shorter 
than ~ 50-100 ms so that they can turn on before collapsing to a black hole. Rotational periods of the hypermassive 
neutron stars are 0.5-1 ms. Thus, if the magnetic fields grow in the dynamical time scale associated with the rotational 
motion via MRI, the amplitude and frequency of gravitational waves may be significantly affected. According to a 
theory of MRI [30], the wavelength of the fastest growing mode is ~ 10(-B/10 12 gauss)(p/10 15 g/cm 3 )^ 1 / 2 (P/l ms) 
cm where B, p, and P denotes a typical magnetic field strength, density, and rotational period, respectively. This 
indicates that a turbulence composed of small eddies (for which the typical scale is much smaller than the stellar 
radius) will set in. Subsequently, it will contribute to a secular angular momentum transport for which the time scale 
is likely to be longer than the growth time scale of MRI ~ a few ms although it is not clear if it is longer than ^100 
ms. On the other hand, if the transport time scale is not as short as ~ 100 ms, other effects associated with magnetic 
fields will not affect the evolution of the hypermassive neutron stars. Indeed, Ref. [12] indicates that the typical time 
scale associated with magnetic braking (winding-up of magnetic field lines) depends on the initial strength of the 
magnetic fields, and it is much longer than the dynamical time scale as ~ 100 (10 12 gauss/B) s. In this case, the 
hypermassive neutron stars can be strong emitters of gravitational waves as indicated in [69]. As is clear from this 
discussion, it is important to clarify the growth time scale of magnetic fields in differentially rotating neutron stars. 
This is also the subject in our subsequent papers [63]. 
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